Molecular junctions in the Coulomb blockade regime: rectification and nesting 
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Quantum transport through single molecules is very sensitive to the strength of the molecule- 
electrode contact. Here, we investigate the behavior of a model molecular junction weakly coupled 
to external electrodes in the case where charging effects do play an important role (Coulomb blockade 
regime). As a minimal model, we consider a molecular junction with two spatially separated donor 
and acceptor sites. Depending on their mutual coupling to the electrodes, the resulting transport 
observables show well defined features such as rectification effects in the I-V characteristics and 
nesting of the stability diagrams. To be able to accomplish these results, we have developed a 
theory which allows to explore the charging regime via the nonequilibrium Green function formalism 
parallel to the widely used master equation technique. Our results, beyond their experimental 
relevance, offer a transparent framework for the systematic and modular inclusion of a richer physical 
phenomenology. 

PACS numbers: 85.65.+h, 73.23.-b, 73.23.Hk, 85.30.Kk 



I. INTRODUCTION 



Single molecule electronicsi'^'^ has been mostly inves- 
tigated in the high temperature and strong contact to 
the electrode regime. The opposite limit of low tempera- 
ture and weakly coupled molecular junctions pose a chal- 
lenge to the currently available experimental techniques. 
Still the possibility to probe the spectroscopy of single 
molecule junctions via a lateral gate could offer new in- 
sights in the peculiar coupling of the electrical and me- 
chanical degrees of freedom at the nanoscale. In order to 
be able to establish the transport mechanisms governing 
such molecular junctions in the Coulomb blockade (CB) 
regime, a technique which could tackle on one hand sin- 
gle electron charging effects and, on the other hand, the 
inclusion of the electron- vibron coupling is of extreme im- 
portance. The nonequilibrium Green function (NEGF) 
formalism has been recently employed to describe trans- 
port observables on the base of a density functional the- 
ory description of the electronic structure^i^i^i^iiiSiiSiiiS 
and model Hamiltonian approaches Jii^iii^ii^ The NEGF 
was applied to describe the influence of the vibron dy- 
namics onto a molecular transistor in the strong coupling 



regim( 



.15.16 



but it is typically substituted with master 



equation approaches when coming to the case of weak 
coupling to the electrodes. Our purpose is to study the 
problem of a two site donor/acceptor molecule in the CB 
regime within the NEGF as a first step to deal with the 
phenomenology of a rigid multilevel island. The nuclear 
dynamics (vibrations) always present in molecular junc- 
tions could be then modularly included in this theory. 
Our method developed in this paper can be calibrated 
on the well-studied double quantum dot problenii2. and 
could be possibly integrated in density functional theory 
based approaches to molecular conductance. 
Here, we apply our theory to the case of a two site ener- 
getically asymmetric molecular junction. In the case of 
serial coupling to the electrodes, this configuration con- 
sists, de facto, in a molecular rectifier (diode) as proposed 



long time ago by Aviram and Ratner— and recently ex- 
perimentally realized We show that the sequential tun- 
neling regime, being a fundamental different regime from 
coherent transport, is compatible to the observed rectifi- 
cation features The serial arrangement of a double site 
correlated molecule between two leads is possibly the sim- 
plest configuration. The most general case (see Fig. [IJ, 
which includes parallel pathways, shows in the sequential 
tunneling regime an interplay of correlated effect and in- 
terference eventually bringing to the phenomenon of a 
nesting of the stability diagrams due to possible different 
charging energies. 

In this paper, we introduce a powerful Ansatz for the 
NEGF which is related both to the equation-of-motion 
(EOM) method and to the Dyson equation approach. 
From the knowledge of the Green function (GF) we then 
calculate the transport observables. Our results are of a 
particular interest in its own at a formal level. In the case 
of a single site junction (SSJ) with Coulomb interaction 
(Anderson impurity model), the conductance properties 
have been successfully studied by means of the EOM 
approach in the cases related to CB^ and the Kondo 
effect 1^ Later the same method was applied to some 




FIG. I: (color) The general configuration of a double site 
junction. The levels ei,2 with charging energies Ui,2 are con- 
nected via t and coupled to the electrodes via the linewidth 
injection rates 7^. 
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two-site modelsiiii22ii^ Multi-level systems were started 
to be considered only recently i^^'^^ Besides, there are 
some difficulties in building the lesser GF in the nonequi- 
librium case (at finite bias voltages) by means of the 
EOM method^i^i^ Here, we develop a self-consistent 
nonequilibrium method for the GF of a single-site junc- 
tion (SSJ) and of a doublc-sitc junctions (DSJ). The re- 
sults of the EOM method could be calibrated with other 
available calculations, such as the master equation ap- 
proach and the non-crossing approximation. 
This paper is organized as follows: after the derivation of 
the SSJ and the DSJ nonequilibrium results for the re- 
tarded and lesser GFs (Sec. II), we do show their effects 
on the transport observables (Sec. III). 

II. SYSTEM AND METHOD 

The goal of this paper is the determination of the trans- 
port observables for a minimal model of a molecular- 
junction in the CB regime, namely a double site cor- 
related impurity Hamiltonian coupled to extended elec- 
trode states. For clarity, we first describe our method in 
the more familiar problem of a single site junction, which 
is the well-known Anderson impurity model. 
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FIG. 2: (color) Tlie stability diagram of a SSJ witli ~ 
2.0 eV, U = 4.0 eV, Fl =Tr = 0.05 eV. (a) The uncorrect 
result obtained by means of the widely used formula in Eq. (|9]) 
for the lesser GF is not symmetric for levels and €cr + U. 
(b) Results obtained by means of our Ansatz in Eq. ([7} shows 
correctly symmetric for levels eo- and + U. 



are the electron self-energies. 



A. Single site case 

The Anderson impurity model is used to describe the 
Coulomb interaction on a single site: 

H = Hd + ^(ifc( + Hau), 
a 

where 

cr 

~ ^ ^ ^k.a'^a.k.a'^a.k.a' 

^aB = ^ {ya,k,a'^,k,a'^a- + ^a,k,<7'^t^a,k,a^ ' 
k,(T 

where d and c are the operators for electrons on the dot 
and on the left {a = L) and the right (q = R) lead, U 
is the Coulomb interaction parameter, is the a level 
of the quantum dot, while e^^ are a level of lead a in fc 
space, a =t, |. With the help of the EOM and a trunca- 
tion approximation, wc can get a closed set of equations 
for the retarded and advanced GFs Gct{t,— 

(uj-e,- K^'')G%^ = S.,r + C/G^'|■■/^ (la) 
{u;-e,-U- S;^/^)Gi')'^/^ - (ns)<5.,., (lb) 

where G^{^ = g'^S"^" = ((nsrfa|4»''/" and 
K^^H = Sf/^ + = y ^ (2) 

a,k '^''^ 



1. Mapping on retarded Green functions: equilibrium case 

There are two typical ways to calculate GFs. The first 
is by means of the Dyson equation and Fcynman dia- 
grams, the second is by means of the EOMi^ For re- 
tarded GFs, from the EOM method, and with the help 
of Eqs. Pa|) and (jlb[) . we can get 

G' = Gl + GIUG^'^'>' 

= g^ + g^se°^g(i>, 

where G"^ is single-particle GF matrix 

" \Gh GlJ ' 

and Gff^r'^ = G^^)r /{ns). Gq describes the single-particle 
spectrum without Coulomb interaction, but including 
the effects from the electrodes. E^^^ = U (rig) is the 
Hartrec-like self-energy of our model. Since there is only 
Coulomb interaction on the site with the levels e^, the 
Fock-like self-energy is vanishing. 

Alternatively, by means of the Dyson equation and the 
second-order truncation approximation, taking Hartree- 
like self-energies S^^. = U{ng) (= E^^^), we can also 
get the retarded GFs as follows;^ 

G^ = Gl + Gl-i:'^G\, (3) 

where G^ = Gq + GqE^Gq is the first-order truncation 
GF. 
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Within the level of the second-order truncation approx- 
imation, we see that there is a map between the EOM 
results and the Dyson results: 



I 



Gq 



I 

Gr /^r 
1 



EOM, (4a) 
Dyson. (4b) 



Eqs. ([4]) prompts a way to include further many-particle 
effects into the Dyson equation. (Eq. (j4b|) . by replacing 
the Dyson-first-order retarded Green function G\ with 
the EOM G^^-''". Then one obtains already the correct 
results to describe CB while keeping the Hartree-like self- 
energy. 



2. Mapping on contour Green functions: nonequilibrium 
case 
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Introducing now the contour GF G, we can get the 

2 

r 



Dyson equation as follows ^9.31,32,33 



G — Go + Go^^jG, 



(5) 



where £ is the self-energy matrix<^ 
According to the approximation for the retarded GF in 
Eq. ([3]), we take the second-order truncation on Eq. (O, 
and then get 

G = Go + GoS'^Gi, 

where Gi = Go + GoS^Go is the first-order contour GF, 

and Go has already included the lead effects. 

Similar to the mapping in Eq. ([?]), we perform an Ansatz 

consisting in substituting the Dyson-first-order G^^^^^ 
with the EOM one 

(^(i)r/a/< consider more many- 
particle correlations, while the EOM self-energy is used 
for the Dyson equation for consistency: 



G 

: 

G 



Go 



Go £^ G\ 
T 

G(i) 



Dyson, 



EOM. 



(6) 



Then, using the Langrcth theorem)^ we get the lesser 
GF, 



G< = G^^ + G^S"''G(i)< + G;fS"'^G(i) 



— Gq 



Gf,f7G(2)< + G;j^C/G(2) 



(7) 



where G, 



•/a/< 



effects, i.e. 



arc GFs for L/ = 0, hut including the lead 



Go< = .9o< + .9oS<Gg + g<^^Gl + g^J:'G<, 
Go" = 9o" + 5S^'S'/"G"■/^ 

with g'f/^^^ the free electron GF. and 



J]r/a/< 



r/a/< 



^r/a/< 



FIG. 3: (color) The stability diagram of serial DSJ with ei,CT = 
e2,<T = -0.15 eV, Ui = U2= 0.3 eY,t 
0.02 eV, 7? = 7r = ,Vi.i.s = 0.005V. 



0.05 eV, jl = 7l 



S< = iE„r,./a(^), and F, = 1(2^^ - S^), Ulu) = 
/(cj — /Iq), / is the equilibrium Fermi function and /i^ 

is the electro-chemical potential in lead a; Sj/^ arc the 
retarded/advanced electron self-energies from Eq. Q and 



G 



(l)r/a/< 



G 



(2)r/a/< 



/{ng-). Performing the same Ansatz 



on the double-particle GF, from Eq. (jlb[) we can get 



G^' 



2)< 



G'(2)rj^(2)<(^(2) 



(8) 



with = S</(ns). 

The lesser GFs in Eq. ([7|) can also be obtained directly 
from the general formula^^ 



G<(w) = G^+ G'oT.'G< + GIY.<G'' + G^^S^G^, 

with the help of the Ansatz in Eq. ([6]). It should be noted 
that Eq. ((T]) is very different from the lesser GF formula, 



G< = g't;<g^ 



(9) 



which is widely used for both first-principle^i2ii2i and 
model Hamiltonian calculations F^i It should be noted 
that the self-energy S< in Eq. ^ contains only con- 
tributions from the electrodes. 

The numerical calculation results of conductance depen- 
dence on the bias and gate voltages by the two differ- 
ent NEGF Eqs. and ^ are shown in Fig. H As we 
can see in the left panel, the adoption of Eq. ([9]) results 
in incorrectly symmetry-breaking in the gate potential. 
This wrong behavior is corrected in the right panel where 
Eq. ([7]) has been used. 
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B. Double site case 

We now investigate the DSJ system with Coulomb in- 
teraction on each dot site. The Hamihonian is expressed 
as follows, 



where 



U, 



3,0 



H. 



E 

k.{T 



with I, j = 1,2 indicate the dot site, t is the constant for 
electron hopping between different sites. 
With the help of the EOM, and by means of the trun- 
cation approximation on the double-particle GFs, we ob- 
tain the closed form for the retarded GFs as follows 



-r) ^ ,.^(2)(C/,t)r ^(C/,t)r 



lUJ 



l,(7-.J,T 

((7,t) 



(10a) 



(10b) 



where the DSJ retarded GFs are defined as G-j.^^'^ 



{{d,Jdl,)Y\ G^j^f' = {{n,^,d.^Jdl,)Y. t means 
'NOT i\ and ^ are the electron self-energy from leads. 
From Eqs. pOap . (jlOb[) and performing the same Ansatz 
as in the case of SSJ, we can obtain the DSJ lesser GFs 
with Coulomb-interaction effects as follows, 



GW*)<(c^) = (1 + G(^'*)^Si-)G(^)< 
•(1 + E^G^'^'*)^) + G^^^^'-^Ef G(^^*'" 



(11) 



with 



fQ t Q 0\ 

f 

t 

vo t oy 



and E< = 0. G^^X is the DSJ lesser GF with the same 
form as Eq. ([7]), but taking 



U = 



^C/i 

;72 

;7i 

t/2y 



, 7^ 

" ~ ' 7I 



where 7^ indicates the line width function of lead a to 
site I, and Ui is the charging energy at site i. G'^^ and 
g.(2)r/a ^Yic GF matrix from Eqs. (fTOa|l and (flObll . 
Here, in order to distinguish different GFs, we introduce 
the subscript '{U,ty for the one with both Coulomb in- 
teraction U and inter-site hopping t, while '(?/)' for the 
one only with Coulomb interaction. 
For our models, the lesser GFs in Eq. © and (fTTjl , 
which are obtained with help of our Ansatz, can also be 
obtained by the EOM NEGF formula in Eq. or in 
Ref. [2^ within the same truncation approximation. 



III. TRANSPORT OBSERVABLES FOR THE 
DOUBLE SITE JUNCTION 

The current can be generally written as^ 

+[/L(u)rL - /H(i,.)r„l(G"'''>' - G""'"")), 

where the lesser GF is given by Eq. pT|) . The differential 
conductance is defined as 



g = 



dJ 



where the bias voltage is defined as Vbias = (/^R ~ A'L)/e. 
A. Serial configuration 



By taking 7^ = 7]^, = 0, wc obtain a serial DSJ, which 
could describe the kind of molecular quantum dot junc- 
tions like the ones in Ref. [Tol. First, at small bias voltages, 
the conductance with the two gate voltages Vg^ and Vg^ 
was calculated, and the relative stability diagram was ob- 
tained as shown in Fig. [3] Because of the double degen- 
eracy (spin-up and spin-down) considered for each site 
and electrons hopping between the dots, there are eight 
resonance-tunnelling regions. This result is consistent 
with the master-equation approachiii 
Further, we studied the nonequilibrium current for large 
bias- voltages (Fig. S]). Because ei^„ and £2,0- are taken 
as asymmetric, for the case without Coulomb interac- 
tion, the I-V curve is asymmetric for ± Vbias, and there 
arc one step and one maximum for the current. The 
step contributes to one peak for the conductance. When 
we introduce the Coulomb interaction to the system, the 
one conductance peak is split into several: two peaks, 
one pseudo-peak and one dip, while the current maxi- 
mum comes to be double split (see Fig.|3|). This process 
can be understood by the help of Fig. [5l At zero bias- 
voltage, £2,0- is occupied and ei^a- is empty, a) By adding 
a bias voltage, the Fermi window is opened. The level 
£2.0- + [/ is first resonant with the edge of the window. 
It will contribute the first peak for conductance, b) By 
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FIG. 4: (color) Current and conductance vs. bias-voltage of 
a DSJ far from equilibrium with parameters ei,o- ~ 0.5 eV, 
£2 <T = -0.5 eV, Ui = U2 = U = 0.2 eV, t = 0.07 eV, jl = 
7i = 0.03 eV, Fg, = -Fg, = V'bia./4 and = -li = 
Vbias/2. The red curve represents the current, while the blue 
the conductance. The inset is the blow-up for the conductance 
peak split. The dash and dot-dash curves are for current and 
conductance with U — 0, respectively. 

opening the window further, the levels e2,cr and ei^o- come 
into the window resulting in the second peak, c) When 
the level ei^a + U comes in, only a pseudo-peak appears. 
This is because there is only a little possibility for elec- 
trons to occupy the level ei^o- under positive bias voltage, 
d) Levels 62,0- -I- U and ei^o- meet, which results in elec- 
tron resonant-tunnelling and leads to the first maximum 
of the current. Then a new level ei^^ + U appears over 
the occupied ei^a for charge effects, e) In Fig. 4(e), the 
meeting of €2,0- and ei^o- results in electron resonant tun- 
nelling. It means that ei^o- will be occupied, which leads 
to the appearance of a new level ci^a + U. Then e2,cr + U 
meets ei ,j -I- U and another resonant tunnelling channel 
is opened for electrons. The two channels result in the 
second current maximum, f) Fig. 4(f) shows that the 
level ei.o- -I- U disappears if the level ei^o- is empty. This 
means that a dip appears in the conductance. 
It should be noted that the characteristics of serial DSJ 
in Fig. [4] have showed some reasonable similarities to ex- 
periments of a single-molecule diodei^^ 

B. Parallel configuration 

If on the other hand, the two sites are symmetrically 
connected to the electrodes, possibly with a small inter- 
dot hopping, but with different charging energies Ui and 
U2, then the stability diagram would show some nesting 
characteristics (Fig. (5]) . 
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FIG. 5: (color) The processes involved in the transport char- 
acteristics in figure m ei = ei,o-, 62 = £2,0-, el = ei,cr -I- U, 
62 = £2,(7 + U. The red line indicates electron resonant- 
tunnelling, a) The first conductance peak, b) The second 
conductance peak, c) The pseudo-peak of conductance, d) 
The first current maximum, and the red line indicates reso- 
nant tunnelling of electrons, e) The second current maximum 
for electron resonant tunnelling, f) The dip of conductance. 



The physics of the thin lines in the figure can be under- 
stood by the help of charging effects. For simplicity, here 
we would ignore the site index i. In the region of large 
positive gate voltage at zero bias voltage, e| and are all 
empty, which means that the two levels are degenerate. 
Therefore adding a bias voltage, first, there will be two 
channels (e^ and e^) opened for current (thick lines). Af- 
ter then, one level (spin-up or spin-down) is occupied, 
while the other obtains a shift for Coulomb interaction: 
eg ^ Eg- + U. Therefore, when the bias voltage is further 
increased to make the Fermi-window boundary meeting 
level eg + U , only one channel is opened for the current, 
which results in the thin lines in Fig. [6l The similar case 
appears in the region of large negative gate voltages. 



IV. CONCLUSIONS 

In this paper, we introduced a powerful Ansatz for the 
lesser Green function, which is consistent with both the 
Dyson-equation approach and the equation-of-motion 
approach. By using this Ansatz together with the 
standard equation-of-motion technique for the retarded 
and advanced Green functions, we obtained the NEGF 
for both the single and the double site junctions in 
the Coulomb blockade regime at finite voltages and 
calculated the transport observables. The method can 
be applied to describe self-consistently transport through 
single molecules with strong Coulomb interaction and 
arbitrary coupling to the leads. 

To test our method, we analyzed the CB stability dia- 
gram for a SSJ and a DSJ. Our results are all consistent 
with the results of experiments and the master-equation 
approach. We showed, that the improved lesser Green 
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FIG. 6: (color) Nested stability diagram of a parallel DSJ 
with parameters ei.o- = —1.8 eV, 62, a = —0.3 eV, Ui — 3.6 eV, 
U2 = 0.6 eV t = 0.001 eV, 7l = 7r = 0.04 eV, 7^ = 71 = 
0.05 eV, = Vg,/2 = VJ2 and Vn = -Vl = Vbia./2. 

function gives better results for weak molecule-to-contact 
couplings, where a comparison with the master equation 
approach is possible. 

For the serial configuration of a DSJ, such as a 
donor/acceptor rectifier, the I-V curves maintain 
diode-like behavior, as it can be already inferred by 
coherent transport calculations 1^ Besides, we predict 
that as a result of the charging effects, one conductance 
peak will be split into three peaks and one dip, and 
one current maximum into two. For a DSJ parallel 
configuration, due to different charging energies on 
the two dot sites, the stability diagrams show peculiar 
nesting characteristics. In both cases, we present the 
results of numerical calculations as well as the simple 
qualitative picture of physical processes. 
We believe, that the results presented here, beyond 
their experimental relevance, might be the transparent 
base for a systematic and modular inclusion of a richer 
physical phenomenology. Work is currently in progress 
to include the electron- vibron interactions to this theory. 
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APPENDIX A: DERIVATION OF THE EOM 
LESSER GREEN FUNCTION 

From the view of perturbation theory, our Hamiltonian 
can be generally written as iJ = Hq + Hi, where Hi is the 
perturbation term to the solved Hq. The contour-ordered 
GF is defined by means of the Schwingcr-Keldysh time 
contour 

{{A{ri);B{r,))f ^ -i(Tc{A(ri)i3(r2)}), (Al) 

where A{ti) and B{t2) are Heisenberg operators, defined 
along the contour C. Taking the time derivative, we ob- 
tain the EOM as, 

+ {{[A{Ti),Hi],B{T2))f. 

Using the free particle solution .g'^(Ti— T2), we can rewrite 
the time-dependent solution as 

{{A{ri);B{r2))f=g''{ri-T2)mTi),B{T2)]±) 
+ 1 9''{n-r'mA{r'),Hi];B{r2))fdr'. 

Now applying the Langreth theorem and transforming in 
the spectral space, we get 

{{A\B))< = g<{cu){[A,B]±)+giu;y{{[A,Hi],B))< 
+ gi^)<{{[A,mlB))Z. (A2) 
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